%{
AUTHOR: Felipe Arteaga

DATE:  2020--
-------------------------------------------------------------------------
PROJECT:
-------------------------------------------------------------------------
DESCRIPTION:

Evalua efecto de SMS enviado el 2017 a postulantes que postularon antes de
que el API fuera activada. El SMS tambien tiene threshold de riesgo que
variaron según lugar

Se aleatorizó entre 4 SMSs:

1) 6 de cada 10 postulaciones como la de JUAN no quedan en los colegios escogidos porque no hay vacantes suficientes. Te recomendamos buscar mas colegios
2) Puede pasar que JUAN no quede en los colegios escogidos porque no tienen vacantes suficientes para todos. Te recomendamos buscar y agregar mas colegios.
3) Muchas familias postulan a los mismos colegios que JUAN y no hay vacantes suficientes para todos. Te recomendamos buscar y agregar mas colegios
4) Te recomendamos agregar mas colegios a la postulacion de JUAN ya que muchas familias estan postulando a los mismos colegios y no hay vacantes suficientes

1) 6 out of 10 applications such as JUAN's are not assigned to any of the chosen schools because there are not enough vacancies. We recommend you search for more schools
2) It may happen that JUAN will not be assigned to any of the chosen schools because they do not have enough vacancies for everyone. We recommend you search and add more schools.
3) Many families apply to the same schools as JUAN and there are not enough vacancies for everyone. We recommend you search and add more schools
4) We recommend adding more schools to JUAN's application since many families are applying to the same schools and there are not enough vacancies.


=========================================================================
%}

clc;clear;close all;fclose('all');
feature('DefaultCharacterSet', 'UTF-8');

pcName=char(java.lang.System.getProperty('user.name'));
if(strcmp(pcName,'felipe'))
    % PC Felipe
    myDir='/Users/felipe/Dropbox/';
    projectDir=[myDir,'git/warnings/'];
    projectDirData=[myDir,'projects/warnings/'];
end


addpath(genpath([myDir,'/myMatlabFunctions/']));

%%
dirPlots=[projectDir,'/paper/figuresCL/playRDs/'];

% warning('Alternative dir plots')
% dirPlots=[projectDir,'/paper/figuresCL/RDs/Aux/'];

dirTable=[projectDir,'/paper/tablesCL/'];
dirData=[projectDirData,'/data/chile/'];

tableBandwidthComparison=false;
plotRDs=false;
titlePlotRD=true;
plotWithZoom=false;
savePlotRDs=false;
bandwidthTable='local'; % 'full','local','localAuto';

switch bandwidthTable
    case 'full'
        bwnote='Full bandwidth was used with 2nd order (local) polynomial fit. ';
    case 'local'
        bwnote='Bandwidth was set to +-0.1 with default local polynomial fit (1st order). ';
    case 'localAuto'
        bwnote='Bandwidth was set to automatic (bwselect with default options) with default polynomial fit (1st order). ';
end

%% Load data:
anho=2017; %0: pooled
if(anho==0)
    anhoStr='';
else
    anhoStr=sprintf('%i',anho);
end


load([dirData,anhoStr,'/inputRD'],'dataRD')


dataRD.rural=dataRD.newMarket<0;

% As asignment at .3 is control, just move (if any) in .3 to epsilon to the
% left:
dataRD.riskSMS(dataRD.riskSMS==.3)=.29999999;
dataRD.riskSMS(dataRD.riskSMS==.5)=.49999999;
dataRD.riskSMS(dataRD.riskSMS==.7)=.69999999;



% Avoid mass points at extremes:
dataRD.restrAll=dataRD.pobSMS==1&dataRD.riskSMS>.01&dataRD.riskSMS<.99;

% No need to be more specific, and put "Predicted placement risk of 1st
% attempt"
dataRD.Properties.VariableDescriptions{'riskSMS'}='Predicted placement risk';
dataRD.Properties.VariableDescriptions{'treatedSMS'}='Received SMS';






%% RD

% Defines subpopulation in which I want to calculate RDs for outcomes
% defined soon.

assert(all(dataRD.cutoffSMS(dataRD.restrAll)>0))



subpobs=cell(1,8);
s=1;
subpobs(s,:)={'Pooled',    dataRD.restrAll,'pooled','addAnyPreSEnd','0.0','.28','.28','riskSMSCentered'};s=s+1;
subpobs(s,:)={'0.3',   dataRD.restrAll&abs(dataRD.cutoffSMS-0.3)<.001,'c.3','','0.3','.28','.68','riskSMS'};s=s+1;
subpobs(s,:)={'0.5',      dataRD.restrAll&abs(dataRD.cutoffSMS-0.5)<.001,'c.5','','0.5','.48','.48','riskSMS'};s=s+1;
subpobs(s,:)={'0.7',    dataRD.restrAll&abs(dataRD.cutoffSMS-0.7)<.001,'c.7','','0.7','.68','.28','riskSMS'};s=s+1;
subpobs(s,:)={'SMS1 - Pooled',    dataRD.restrAll&not(ismember(dataRD.tipoSMS_A,[2 3 4])),'SMS1_pooled','','0.0','.28','.28','riskSMSCentered'};s=s+1;
subpobs(s,:)={'SMS2 - Pooled',    dataRD.restrAll&not(ismember(dataRD.tipoSMS_A,[1 3 4])),'SMS2_pooled','','0.0','.28','.28','riskSMSCentered'};s=s+1;
subpobs(s,:)={'SMS3 - Pooled',    dataRD.restrAll&not(ismember(dataRD.tipoSMS_A,[1 2 4])),'SMS3_pooled','','0.0','.28','.28','riskSMSCentered'};s=s+1;
subpobs(s,:)={'SMS4 - Pooled',    dataRD.restrAll&not(ismember(dataRD.tipoSMS_A,[1 2 3])),'SMS4_pooled','','0.0','.28','.28','riskSMSCentered'};s=s+1;
subpobs(s,:)={'Pooled R12 conge',    dataRD.restrAll&(dataRD.newMarket==12001|(dataRD.newMarket==12002&dataRD.grade<9))&dataRD.lengthPreS<5,'congR12','addAnyPreSEnd','0.0','.28','.28','riskSMSCentered'};s=s+1;
%subpobs(s,:)={'R12 (0.3)',    dataRD.restrAll&abs(dataRD.cutoffSMS-0.3)<.001&(dataRD.newMarket==12001|dataRD.newMarket==12002),'r12','','0.3','.28','.68','riskSMS'};s=s+1;
%posAll=1;




% Definition of outcomes for RDs. 3rd row defines the beginning of a panel
% in the table
% 4th row


% New: enrolled conditional on placed (short name because of stata)
dataRD.enrolledInAssignedCOA=double(dataRD.enrolledInAssigned);
dataRD.enrolledInAssignedCOA(dataRD.assignedToPref==0)=nan;
dataRD.placedInAnyPreferenceNWIE=dataRD.placedInAnyPreferenceNoWaitlistIniEnd;
dataRD.placedInAnyPreferenceNoRiskIE=dataRD.placedInAnyPreferenceNoRiskIniEnd;

%notIn2020=[notIn2020,'enrolledInAssignedCOA'];

% 5th column defines if IV is calculated or not.
outcomes={'changeAnyIniEnd','Any modification','cha','',0;...
    'addAnyPreSEnd','Add any','aa','',0;...
    %'deltaProbPlacedNoRiskIniEnd','$\Delta$ prob. placed to zero risk','dpnr','',1;...
    %'deltaProbPlacedNoWaitlistsIniEnd','$\Delta$ prob. placed to uncongested','dpw','',1;...
         'schoolsAddedPreSEnd','Schools Added','sa','',1;...
    %     'riskLastDayIni','True Risk initial attempt','ri','';...
    %     'riskLastDayEnd','True Risk final attempt','rf','';...
         'DRiskExpostPreSEnd','$\Delta$ Risk','drEP','',1;...
    %     'riskExpostIni','True Risk initial attempt','riEP','';...
    %     'riskExpostEnd','True Risk final attempt','rfEF','';...
    %     'DRiskExpostIniEnd','$\Delta$ Risk','drEP','',1;...
    %     'numberOfNewUncongestedIniEnd','Uncongested schools Added','saUncong','';...
    %     'addAsFirstIniEnd','Add as first','af','',1;...
    %     'addInBetweenIniEnd','Add to middle','ab','',1;...
    %     'addAsLastIniEnd','Add as last','al','',1;...
    %     'deleteAnyIniEnd','Drop any','da','',1;...
    %     'changeOrigOrderIniEnd','Re-order','ro','',1;...
    %     'riskLastDayEnd','Risk final portfolio','rf','';...
         'placedInAnyPreferencePreSEnd','Placed to preference','ap','',1;...
    %     'placedInAddedIniEnd','Placed to added preference','apad','';...
    %     'placedInOldIniEnd','Placed to old preference','apold','';...
    %     'enrolledInAssigned','Enrolled in placed','ea','',1;...
    %     'enrolledInAssignedCOA','Enrolled in placed$|$placed','eaadcp','',1;...
    %     'declineOffer','Decline offer','do','',1;...
    %     'enrolledInAddedIniEnd','Enrolled in added','eaad','';...
    %     'addAnyUncongestedIniEnd','Add any uncongested','aaUncong','D. Congestion-related outcomes',1;...
    %     'placedInAnyPreferenceNWIE','Placed to uncong.','apUncong','',1;...
    %     'placedInAddedNoWaitlistIniEnd','Placed to uncong. added pref.','apadUncong','',1;...
    %     'enrolledInAddedNoWaitlistIniEnd','Enrolled in uncong. added','eaadUncong','',1;...
    %     'addAnyNoRiskIniEnd','Add any zero risk','aaNR','D2. Risk-related outcomes',1;...
    %     'placedInAnyPreferenceNoRiskIE','Placed to zero risk','apNR','',1;...
    %     'placedInAddedNoRiskIniEnd','Placed to zero risk added pref.','apadNR','',1;...
    %     'enrolledInAddedNoRiskIniEnd','Enrolled in zero risk added','eaadNR','',1;...
    %     'acceptOffer','Accept offer','ao','';...
    %     'preferenciaAsign_1','Order assigned','oa','';...
    };



if(anho==2020)
    outcomes=outcomes(not(ismember(outcomes(:,1),notIn2020)),:);
end



dataRD.Properties.VariableDescriptions{'riskSMS'}='Predicted placement risk';
for o=1:size(outcomes,1)
    
    dataRD.Properties.VariableDescriptions{outcomes{o,1}}=outcomes{o,2};
end


assert(allunique(outcomes(:,1)))
assert(allunique(outcomes(:,2)))
assert(allunique(outcomes(:,3)))

dataRD.riskSMSCentered=dataRD.riskSMS-dataRD.cutoffSMS;
data=dataRD(:,[{'riskSMS','cutoffSMS','riskSMSCentered','tipoSMS_A'},outcomes(:,1)']);
inAnyPob=false(height(data),1);


% Generate stata commands
commands=cell(size(subpobs,1)*size(outcomes,1),2);
counter=1;


for p=1:size(subpobs,1)
    

    cutoffStr=subpobs{p,5};
    lbwStr=subpobs{p,6};
    rbwStr=subpobs{p,7};
     runVar=subpobs{p,8};
    
    data.(sprintf('subpop%i',p))=subpobs{p,2};
    inAnyPob=inAnyPob|subpobs{p,2};
    
    for o=1:size(outcomes,1)
        
        % Avoid computing outcomes that are not available:
        notIn2020={};
        noCalcular=(strcmp(subpobs{p,3},'2020')|strcmp(subpobs{p,3},'2020comp'))&ismember(outcomes{o,1},notIn2020);
        
        if(not(noCalcular))
            
            
            % Full bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'full')||plotRDs)
                commands(counter,:)={sprintf('%s_%i_full',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, c(%s) p(2) h(%s %s)',outcomes{o,1},runVar,p,cutoffStr,lbwStr,rbwStr)};counter=counter+1;
            end
            % Full loccal fixed bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'local'))
                commands(counter,:)={sprintf('%s_%i_local',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, c(%s)  h(.1 .1)',outcomes{o,1},runVar,p,cutoffStr)};counter=counter+1;
            end
            % Full local auto bandwidht
            if(tableBandwidthComparison||strcmp(bandwidthTable,'localAuto'))
                commands(counter,:)={sprintf('%s_%i_localAuto',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, c(%s) ',outcomes{o,1},runVar,p,cutoffStr)};counter=counter+1;
            end
            
            % If with IV:
            if(outcomes{o,5}==1&&~isempty(subpobs{p,4}))
                
                endogenousVar=subpobs{p,4};
                
                % Full bandwidth
                if(strcmp(bandwidthTable,'full'))
                    commands(counter,:)={sprintf('%s_%i_full_iv',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, fuzzy(%s) c(%s) p(2) h(%s %s)',outcomes{o,1},runVar,p,endogenousVar,cutoffStr,lbwStr,rbwStr)};counter=counter+1;
                end
                % Full loccal fixed bandwidth
                if(strcmp(bandwidthTable,'local'))
                    commands(counter,:)={sprintf('%s_%i_local_iv',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, fuzzy(%s) c(%s)  h(.1 .1)',outcomes{o,1},runVar,p,endogenousVar,cutoffStr)};counter=counter+1;
                end
                % Full local auto bandwidht
                if(strcmp(bandwidthTable,'localAuto'))
                    commands(counter,:)={sprintf('%s_%i_localAuto_iv',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, fuzzy(%s) c(%s) ',outcomes{o,1},runVar,p,endogenousVar,cutoffStr)};counter=counter+1;
                end
            end
        end
    end
end
commands=commands(not(cellfun(@isempty,commands(:,1),'UniformOutput',true)),:);


if(false)
    % Run specific rd
    dataRD.auxRest=dataRD.restrAll&dataRD.mandatoryToAdd==0;
    a=stataCommand(sprintf('rdrobust lengthEnd riskSMS if auxRest==1, c(%s) p(1) h(.1 .1)',cutoffStr),dataRD)
    stataExpress(sprintf('rdrobust  DRiskExpostIniEnd riskSMS if auxRest==1, fuzzy(addAnyIniEnd) c(%s) p(1) h(.1 .1)',cutoffStr),data)
end

% Run Stata:
data=data(inAnyPob,:);
res=stataCommand(commands,data);


%% Make plots and table
close all

cantPobs=(size(subpobs,1));
cantRegs=(size(outcomes,1));

matrix_nl=nan(cantRegs,cantPobs);
matrix_nr=nan(cantRegs,cantPobs);
matrix_b=nan(cantRegs,cantPobs);
matrix_se=nan(cantRegs,cantPobs);

matrix_b_iv=nan(cantRegs,cantPobs);
matrix_se_iv=nan(cantRegs,cantPobs);
matrix_nl_iv=nan(cantRegs,cantPobs);
matrix_nr_iv=nan(cantRegs,cantPobs);

matrix_b_l=nan(cantRegs,cantPobs);
matrix_se_l=nan(cantRegs,cantPobs);

matrix_b_r=nan(cantRegs,cantPobs);
matrix_se_r=nan(cantRegs,cantPobs);

matrix_biasCorr_ci=nan(cantRegs,cantPobs,2);



regNames=fieldnames(res);

for p=1:cantPobs
    for r=1:cantRegs
        
        regName_i=sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable);
        if(ismember(regName_i,regNames))
            res_i=res.(regName_i);
            
            % Beta
            matrix_b(r,p)=res_i.tau_cl;
            % Sd Beta
            matrix_se(r,p)=res_i.se_tau_cl;
            
            % Beta left
            matrix_b_l(r,p)=res_i.beta_p_l(1);
            % Sd Beta left
            matrix_se_l(r,p)=sqrt(res_i.V_cl_l(1));
            
            % Beta right
            matrix_b_r(r,p)=res_i.beta_p_r(1);
            % Sd Beta right
            matrix_se_r(r,p)=sqrt(res_i.V_cl_r(1));
            
            % Bias corrected Confidence interval (alpha=5%)
            matrix_biasCorr_ci(r,p,1)=res_i.tau_bc-norminv(.975)*res_i.se_tau_rb;
            matrix_biasCorr_ci(r,p,2)=res_i.tau_bc+norminv(.975)*res_i.se_tau_rb;
            
            % N
            matrix_nl(r,p)=res_i.N_h_l;
            matrix_nr(r,p)=res_i.N_h_r;
            
            % If with IV:
            if(outcomes{r,5}==1&&~isempty(subpobs{p,4}))
                regName_i_iv=sprintf('%s_%i_%s_iv',outcomes{r,3},p,bandwidthTable);
                
                res_i_iv=res.(regName_i_iv);
                % Beta
                matrix_b_iv(r,p)=res_i_iv.tau_cl;
                % Sd Beta
                matrix_se_iv(r,p)=res_i_iv.se_tau_cl;
                
                % N
                matrix_nl_iv(r,p)=res_i_iv.N_h_l;
                matrix_nr_iv(r,p)=res_i_iv.N_h_r;
                
            end
            
            if(plotRDs)
                
                
                res_i=res.(sprintf('%s_%i_%s',outcomes{r,3},p,'full'));
                res_i_table=res.(sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable));
                
                plotsub=subpobs{p,2};
                
                figure
                plotRD(res_i,dataRD,plotsub,'otherPointEstimate',res_i_table);
                
                if(titlePlotRD)
                    title(sprintf('RD %s - %s',outcomes{r,2},subpobs{p,1}))
                    
                end
                
                if(savePlotRDs)
                    easyExport([dirPlots,sprintf(sprintf('%s_SMS_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
                
                if(plotWithZoom)
                    % Plot with zoom
                    figure
                    plotRD(res_i,dataRD,plotsub,'newxlim',[.1,.5],'nb',100,'otherPointEstimate',res_i_table);
                    easyExport([dirPlots,sprintf(sprintf('%s_SMS_%s_%s_withZoom.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
          
            else
                fprintf('Ojo: reg %s no existe para subpob %s\n',regName_i,subpobs{p,1});
            end
        end
    end
    
    %vector_Nl(1,p)=sum(dataRD.riskSMS<.3&subpobs{p,2});
    %vector_Nr(1,p)=sum(dataRD.riskSMS>.3&subpobs{p,2});
    
end





%% RD main table
matTable=nan(size(matrix_b).*[1 2]+[2 0]);
matTable_sd=nan(size(matrix_b).*[1 2]+[2 0]);

matTable(:,1:2:end)=[matrix_b;...
    matrix_nl(1,:)...
    ;matrix_nr(1,:)];

matTable(:,2:2:end)=[matrix_b_iv;...
    matrix_nl_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)...
    ;matrix_nr_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)];


matTable_sd(:,1:2:end)=[matrix_se;nan(2,size(matrix_se,2))];
matTable_sd(:,2:2:end)=[matrix_se_iv;nan(2,size(matrix_se_iv,2))];

withInfo=not(all(isnan(matTable),1));


header=repmat({''},2,size(matTable,2));
header(1,1:2:end)=subpobs(:,1)';
header(1,2:2:end)=subpobs(:,1)';
header(2,2:2:end)={'IV'};


matTable=matTable(:,withInfo);
matTable_sd=matTable_sd(:,withInfo);
header=header(:,withInfo);

header=[{'Risky cutoff';''},header];

opt=struct;
opt.header=header;
opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
opt.primeracolumna=[outcomes(:,2);{'NL';'NR'}];
opt.addColumnNumber=true;

opt.filaFantasma=size(matTable,1)-2;
paneles=find(not(ismissing(outcomes(:,end-1))));
opt.panel=[num2cell(paneles-1),outcomes(paneles,end-1)];
%opt.alignmentFirstCol={'L{3cm}'};
opt.adjust=true;

opt.file=sprintf('%s_tableSMS%s',dirTable,anhoStr);

opt.title='RD estimates of SMS effects - 2017';
opt.label='tab:RD-SMS2017';

opt.note=['Local linear RD estimates of SMS effects from warning text message. Computed using triangular kernel with bandwidth 0.1. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}. ',...
     'We report estimates in the pooled sample and for each different risky cutoff definition. IV (column 2) shows the instrumental variable specifications, where the endogenous regressor is the add any school indicator. ',...
     ''];

tabla=cell2latex(mat2cellstr(matTable,'precisionDecimal','%.3f','revisarFilas',true),'opts',opt);
compileLatex(tabla);

barmean(data(data.tipoSMS_A>0&data.subpop1,:),'tipoSMS_A','addAnyPreSEnd')
title(sprintf('Only %s',subpobs{1,1}))
figure
barmean(data(data.tipoSMS_A>0&data.subpop9,:),'tipoSMS_A','addAnyPreSEnd')
title(sprintf('Only %s',subpobs{9,1}))
